EfRcient wavefunction propagation by minimizing accumulated action 
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This paper presents a new technique to calculate the evolution of a quantum wavefunction in a 
chosen spatial basis by minimizing the accumulated action. Introduction of a finite temporal basis 
reduces the problem to a set of linear equations, while an appropriate choice of temporal basis 
set offers improved convergence relative to methods based on matrix exponentiation for a class of 
physically relevant problems. 



Calculating the time evolution of a quantum wavefunc- 
tion is a longstanding problem in computational physics. 
The fundamental nature of the problem makes it resistant 
to simplification. Problems of physical interest, such as 
the interaction of a molecule with a strong laser field, may 
lack symmetry or involve time dependent, nonperturba- 
tive fields. Problems involving multiple dimensions or 
multiple interacting particles may quickly grow so large 
as to be unmanageable with all but the largest compu- 
tational resources [l|-l3|, a problem which is not easily 
outstripped by increases in computational power. An 
ideal propagator, then, must serve two masters - it must 
treat the physical side of the problem accurately, and the 
computational side of the problem efficiently. 

The most common approach to the problem does not 
treat the evolving wavefunction directly. Instead, the 
wavefunction 'ip[x,t) is expanded in some chosen basis set 
't>{x,t) = I]iCi(t)xj(x), where Ci(to) = (XilV'(^o))- Here 
ip{x,t) is the true wavefunction and (f>{x,t) is its repre- 
sentation in the chosen basis. After this expansion has 
been made, the propagation scheme may operate only on 
the wavefunction's representation, rather than the wave- 
function itself. Within this general framework, there 
has been a great profusion of methods for calculating 
the evolution of the coefficients Ci{t). Popular methods 
include Crank-Nicholson^^], second order differencing ^5,], 
split operator 0, short iterative LanczosQ, and Cheby- 
shev propagation [8|, as well as many others [9[. 

As may be inferred from the large number of com- 
peting methods, the practical question of which method 
works best is very difficult to answer, and usually requires 
problem-specific information. Ironically, the time depen- 
dent Schrodinger equation (TDSE), while challenging to 
solve, is not difficult to satisfy at a particular time. All 
of the above methods satisfy the TDSE exactly or ap- 
proximately at the initial time in the propagation inter- 
val. Indeed, entire families of propagators may satisfy 
the TDSE at the initial point: for the class of propaga- 
tors ip{x, t){l - iHaAt) = ip{x, t + Ai)(l + iH{l - a) At), 
< At, the TDSE at time t is satisfied to first order 
for any value of a. a = 1 yields the forwards Euler 
method, a = the backwards Euler, a = 1/2 the Crank- 
Nicholson. Although such methods may show sharp dif- 
ferences in suitability to particular problems, the TDSE 



alone gives little guidance. A full diagonalization of the 
Hamiltonian matrix would satisfy the TDSE at all times; 
however, such a diagonalization would be prohibitively 
expensive for a large problem, and would not exist for a 
problem with a time dependent Hamiltonian. Compari- 
son of different propagators has often involved numerical 
testing on simple problems 10|, [ll | or algorithmic scaling 
arguments [13 ■ 

This paper addresses the problem of wavefunction 
propagation from the physical perspective of minimiz- 
ing the action accumulated over the chosen time interval. 
Minimizing this action is shown to be equivalent to mini- 
mizing the time integrated error of propagation. Because 
the action is calculated over the entire time step rather 
than at a single point, the constraint that it be minimized 
is more strict than the TDSE, allowing the construction 
of a unique, variationally optimum propagator for a par- 
ticular order in time. 



ERRORS OF PROPAGATION AND 
REPRESENTATION 

A central difficulty of any numerical propagation 
scheme is that, although the propagator seeks to model 
the evolution of an ideal wavefunction ip{x,t), it has ac- 
cess only to the representation of the wavefunction in 
some chosen basis, 4>{x,t). The error of the representa- 
tion is given by S{x, t) = ip{x, t) — 0(a;, t). 

As an alternative to direct exponentiation of the 
Hamiltonian, a propagator may be constructed by mini- 
mizing the integral of the error over the time step 



Global error = 
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Writing the error as a two term Taylor series, 
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and recalling that i-^'tp{x,t) — HTp{x,t) for the true 
wavefunction, the second term in equation [2] can be writ- 



ten as 
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From equations [2] and [3l it is apparent that the global 
error may arise either from imperfectly representing the 
wavefunction in a particular basis (terms containing 
S{x,t)) or from imperfectly describing the evolution of 
the wavefunction in that basis (terms containing 0(x, t)). 
The representation error may be minimized by an appro- 
priate choice of basis; here the focus is on minimizing the 
error of propagation. 

The quantity {4'{x, t)\ (i^ — H) \(j){x^ t)) found in equa- 
tion |3] is the Lagrangian density, and its integral over time 
gives the action accumulated by the wavefunction in a 
particular interval. However, unlike the true Lagrangian 
density, here the action is calculated with respect to the 
representation of the wavefunction, rather than the wave- 
function itself. The distinction is significant. For the 
true wavefunction, minimizing the action is equivalent 
to setting i-^^{x,t) — Hijj{x,t) = for all x and t. For 
the action minimizing representation of the wavefunction, 
I {(t){x,t)\ (id/dt — H) \(t){x,t)) I is dependent on the choice 
of spatial and temporal basis functions and is not guar- 
anteed to be zero. 

For a finite basis of spatial Xi{^) ^-nd temporal r„(t) 
basis functions, a time dependent representation of the 
wavefunction can be written as 
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In this basis, the global error of Eq. [3]becomes a function 
of the coefficients Cm and the matrix representations of 
the quantum mechanical operators. Writing the Hamilto- 
nian as the sum of time independent and time dependent 
operators 



H = Ha{x) + V{x,t) 
and defining the matrices 
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Q„,„ = / dtT* {t)T„,{t), 

the change in action resulting from Ci„ 
given by 



(10) 
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and the condition to minimize the accumulated action is 
that either 6*^ — (for the initial conditions) or 



^i^jm — / ^ ^inl^^ij^nr 



Hij Un 



^ijn7n\^jfn — ^ 



(12) 
for all j,m. In these equations, i appearing as a subscript 
is treated as an index, while i multiplying OijQnm is the 
square root of negative one. 

Equation 1121 is the main result of this paper. In order 
to construct a least action propagator, it is necessary 
only to choose an appropriate temporal basis. While in 
principle this analysis applies equally well to any choice 
of basis, an obvious choice is for T„(i) to be a set of 
linearly independent low-order polynomials in t. 

Lagrange interpolating polynomials provide a conve- 
nient set of temporal basis functions. For an evenly 
spaced grid t„i — t + At * m/n for m ~ 0,n, the in- 
terpolating polynomials are given by 



r™(i) = n 



k—OjU^k^m 



t — tk 



in 



^k 



(13) 



This yields a basis set of n linearly independent n-order 
polynomials in t, with the property that <j){x,tn) — 
J2i ^inXi{x) ■ One advantage of this choice of basis is 
that for small propagation times. Cm will have compa- 
rable amplitudes for all n, making the associated linear 
system easier to solve with high accuracy. 

Having chosen a temporal basis, coefficients Cm which 
satisfy Equation [T^] as well as the initial condition Cjq = 
(XilV'(io)) can be found using Lagrange multipliers. If S 
is the action accumulated in the time interval, let S" = 
S + Y.^^^f*, where /, = Qo - {x^W{to)). The least 
action coefficients are found by minimizing S" with the 
constraint that fi = for all i. This yields a system of 
linear equations 



/ ^ Cin [I'CijQ^ 



for all j,m, and 
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for all i. The Lagrange multipliers Xj calculated in this 
procedure are not needed by the propagator and can be 
discarded after solving the linear system. 



COMPARISON WITH THE LANCZOS 
PROPAGATOR 



The least action propagator derived in the previous 
section is the unique, variationally optimum propagator 
for a particular order in time. As such, it represents 
a formal improvement over all propagators approximat- 
ing the wavefunction as a low order polynomial in time 
- forward and backwards Euler, Crank Nicholson, sec- 
ond order differencing, etc. However, it is less clear 
how this formal improvement translates to a practical 
benefit, or how the least action propagator compares to 
methods which attempt to diagonalize the Hamiltonian 
in a Krylov subspace, such as the popular short iterative 
Lanczos method [7|. 

The Lanczos method works by repeatedly multiply- 
ing the initial wavefunction by the Hamiltonian ma- 
trix to create a Krylov space of limited dimension in 
which the matrix exponential e~'^* can be calculated 
exactly. It is considered to be both efficient and quickly 
converging lOj. Existing variational propagators have fo- 
cused on the evolution of the wavefunction in the Krylov 
subspace, yielding convergence properties similar to the 
Lanczos propagator |13l Il4| . The Chebyshev propagator, 
which also uses repeated multiplication by the Hamilto- 
nian matrix to construct a Krylov space, converges sim- 
ilarly to the Lanczos methodjl5|. 

In the limit that the Krylov space has dimension equal 
to the full Hamiltonian, the Lanczos method is equivalent 



where da^p = Cq, {gplfa), then the error is given by 



to diagonalizing the Hamiltonian, yielding L 



' dt 



^ 



Hip = at all times. This solution is the global minimum 
action solution, and cannot be improved upon. However, 
for most applications, the Krylov subspace is chosen to 
have a much smaller order - typically in the range 1-10. 

Because the Krylov subspace is constructed through 
repeatedly multiplying an initial wavefunction by the 
Hamiltonian matrix, later Krylov basis vectors will tend 
to some overlap with those eigenvectors of H with the 
largest eigenvalues. For problems with Coulomb singu- 
larities or fine spatial bases, it is not uncommon for these 
largest eigenvalues to be artifacts of the choice of ba- 
sis, having no counterpart in the system being described. 
However, they may nonetheless serve to limit the stepsize 
which may be taken with high accuracy. 

If an ideal wavefunction (ie, without reference to a ba- 
sis) "0(2;, t) can be expanded in terms of energy eigenfunc- 
tions over a short time interval 

V^(a:,i)=^c„/„(a;)e-'^-(*-*«) (16) 



and the evolution of the wavefunction's representation in 
the Krylov subspace is given by 
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En — Efl 



-(i-to)). 



(18) 



(19) 



If the Krylov subspace is now partitioned into a "good" 
subspace with Ep < EcutoS and a "bad" subspace with 
Ep > EcutoS, the error can be estimated by setting E^ — 
-E^ = in the good subspace and Ep — Ea — Eh, where 
Eh reflects the largest eigenvalues of H, yielding 

{5\5)^ Y, Mo,/^p4sin2(ii?^(t-<o)). (20) 

Q,/3,B^>_Ecutotf 

The error of the Lanczos method is thus minimized either 
when the projection into the bad subspace is small or 
when En^t « 1. 

In contrast to the Lanczos method, the error of the 
least action propagator is bounded by the error of the 
Taylor series of the true wavefunction. Thus 



{6\S) < J2 K 



N, 



2u-i_E„(t-to) 



E {-iEaY , VI |2 
T^ (*-^o) 



n=0 



and the condition for the error to remain small is simply 
that _EcutoffAi < 1. For a problem with i^cutoff << Eh, 
the least action propagator offers the possibility of much 
larger stepsizes at high accuracy than the Lanczos prop- 
agator. 

The two propagators were tested numerically using a 
1 dimensional Coulomb potential \/x for x ranging from 
to 10. The region was separated into 100 finite el- 
ement regions, with two quadratic finite elements per 
region. The wavefunction was restricted to have zero 
value at both endpoints. The largest eigenvalue of the 
resulting Hamiltonian matrix was 499 Hartree. The ini- 
tial wavefunction was chosen to be a Gaussian of unit 
width, centered at x=2. The choice of initial wavefunc- 
tion and potential were made to ensure that the wave- 
function would be far from equilibrium and have a strong 
interaction with the Coulomb potential, as for an electron 
wavepacket scattering from a positive ion. 

The accuracy of the Lanczos- and least action propa- 
gator was calculated by propagating the initial wavefunc- 
tion a single timestep and comparing the resulting wave- 
function with the "true" wavefunction found by directly 
diagonalizing the full Hamiltonian. yielding an error 



\err\ = {5{t + 5t)\5{t + 5t)) . 



(22) 



For the least action propagator, the order of propagation 
is one less than the degree of the polynomial basis func- 
tions. For the Lanczos propagator, the order is given by 



- 


• 






•^^^^P^ 








W, 


^' 


- 




' J^^^ 


1^^^ 


JJy 














/. 


■fyw- 


■^ 




1st order 

3rd order 

5th order 




/ 


/A 






7th order 

— 9th order 











1 1 1 1 1 1 






"~~^\ ~\ 


-0.02 


- 


\ \ 


-0.04 


- 




\ ^ 






1st order 






3rd order 


\ / 


-0.06 


- 


5 th order 

7th order 


\ - 






9th order 


1,1, 


nns 




1 



-1 



FIG. 1. (Color online) Point error (J(a;,i + Ai)|5(x, t + At)) 
of Crank Nicholson (dashed line) , Lanczos (dotted lines) and 
least action (solid lines) propagators, as a function of step 
size. For large stepsizes, the least action propagator is many 
times more accurate than the Lanczos propagator of the same 
order. 



the dimension of the Krylov subspace, starting with for 
the initial wavefunction. 

The error as a function of order and stepsize for the 
two methods is shown in Figure [T] Also shown in the fig- 
ure is the error vs time for the popular Crank Nicholson 
first order propagator [J] . As opposed to the global error 
which was used in the derivation of the least action prop- 
agator, these figures show the point error after a single 
propagation step. 

These results show that the least action propagator of- 
fers the potential for large timesteps to be taken with 
high accuracy, with the greatest advantage coming from 
propagation at high order. At first order, the least action 
propagator gives approximately the same point error as 
the Crank Nicholson method, while higher orders rapidly 
decrease the error for a particular timestep, or alterna- 
tively increase the size of the timestep which can be taken 
for a particular desired error. As the order increases, the 
error begins to saturate as different order propagators 
converge on the same result. That this saturation does 
not result in zero error may result from numerical error 
in the diagonalization of the Hamiltonian or the linear 
solver. 

For all orders tested, the least action propagator was 
many times more accurate than the Lanczos propagator 
for large stepsizes. For very small stepsizes, the error of 
both methods was comparable, with the Lanczos method 
more accurate. Both methods became much more accu- 
rate for stepsizes of less than 10~^, which is interpreted 
to mean that the initial wavefunction had some projec- 
tion onto very high energy eigenstates; ie, the sample 
problem did not have -Ecutoff << E^. 

One weakness of the least action propagators which 



FIG. 2. (Color online) Exponential growth rate 

(logjQ |norm|)/di vs propagation time for different or- 
ders of the least action propagator. Higher orders show a 
growth rate closer to zero. 



arises from the choice of polynomial basis functions is 
that the norm of the propagated wavefunction is not 
required to be a constant as a function of time. Fig- 
ure [2] shows the rate of growth/decay of the norm {(\>\(\>) 
for different orders of propagation as a function of step 
size. Here the largest deviation from zero growth rate 
is found for the combination of low order and large step- 
size. Higher order propagators show growth rates close to 
zero. For problems requiring repeated use of the prop- 
agator over many timesteps, it is thus likely that the 
propagated wavefunction will need to be renormalized 
periodically. Because the norm is still very close to 1, 
such renormalization has a minimal effect on the point 
error shown in Figure [TJ 

From these figures, it is apparent that the least ac- 
tion propagator works best at high orders, which offer 
the combination of large time steps, high accuracy and 
low rates of growth or decay of the norm. However, high 
order also increases the size of the linear system which 
must be solved at each step. For a basis consisting of Ux 
spatial and nj temporal basis functions, the least action 
propagator requires rig — n^int-V^^ variables to be solved 
for. While specific implementations are beyond the scope 
of this paper, the question of how best to solve this linear 
system will play a crucial role in applying the least action 
propagator to real world problems. To this end, a few fea- 
tures of the least action linear system are worth noting. 
First among these is the highly separable nature of the 
least action linear system defined in equations [Til and [T5] 
In equation [Ml the variation of the action is given as the 

Of 



sum of three matrices, OijQnm, HijUmm and Vi 



%jnm ■ 



these, the first two are separated into the product of spa- 
tial and temporal matrices. Because of this, these two 
matrices inherit the sparsity and/or banded structure of 
the underlying spatial matrices. The case is similar for 



the nonseparable Vijnm '■ the integral over time and space 
is nonzero only if the integral over space is nonzero. Be- 
cause of this, basis sets such as finite elements which are 
chosen for the structure of their Hamiltonian and overlap 
matrices will retain these advantages in the least action 
equation. For a banded problem such as the ID finite 
element problem treated in this paper, the bandwidth of 
the linear system increases linearly with the order of the 
propagator, giving an overall n^ scaling with the order. 

For very large problems which lack such a simple struc- 
ture, it is likely that solution of the least action linear sys- 
tem will require use of an iterative solver, such as those 
available in the PETSc [il| or Trilinos [i3l libraries. Such 
solvers, require calculating a matrix vector product at 
every iteration. Here the advantages of the least action 
equation's separable form are very apparent, particularly 
in the case of a static Hamiltonian. A single matrix vector 
product of the linear system defined in equations [14] and 
1151 requires 1 (very expensive) matrix- vector multiplica- 
tion by the unseparated matrix Vijnrm 2nf independent 
(expensive) multiplications of the form Dim = MijCjm, 
where M = O or H, followed by 2nx (cheap) indepen- 
dent multiplications of the form Fin = NnmDim, where 
N = Q 01 U. Thus, although the linear system which 
must be solved is very large, it is well suited to itera- 
tive solution. The overall scaling of one iteration with 
respect to propagator order will be limited by the slow- 
est of these steps, which may depend on specifics of the 
data structure and architecture of the system used. 

This paper has addressed the problem of propagating 
a wavefunction in time by minimizing he accumulated 
action. For a particular choice of spatial and tempo- 
ral basis functions, the problem is reduced to solution 
of a (potentially very large) system of linear equations. 
This linear system inherits the sparsity and/or banded 
structure of the spatial Hamiltonian and overlap matri- 
ces, while its separable structure makes it amenable to 
solution by iterative solvers. The resulting propagator 
was shown to have improved convergence relative to the 
commonly used short iterative Lanczos propagator, giv- 



ing the potential for larger stepsizes at high accuracy. 

The derivation of the action as a measure of propaga- 
tion error is a powerful result which offers many oppor- 
tunities for the systematic improvement of propagation 
schemes. By monitoring the spacetime volumes where 
the most action is accumulated, spatial or temporal bases 
could be selectively refined to increase the total accuracy 
of a propagation step at minimal additional computa- 
tional effort. For this reason, the full power of the action 
minimization principle may be yet to be unlocked. 
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